Method of determining earth elastic parameters in anisotropic media

ABSTRACT

A method of determining the elastic moduli of anisotropic earth strata uses acoustic and density data without making simplifying approximations of the Christoffel relation or the assumption of weak anisotropy. To this end, the relation is converted algebraically into a linear system and solved analytically to generate a solution for estimating elastic moduli. Preferably, a set of qP phase points obtained from walkaway VSP data, with multiple vertical angles, plus at least one vertical SV point, are used to derive the elastic moduli A 11 , A 13 , A 33 , A 55 . In the case of FTIV or other orthorhombic media, the method may also be applied to obtain additional elastic moduli.

FIELD OF THE INVENTION

This invention relates to a method of determining elastic parameters ofearth strata, and more specifically to an improved method of determiningelastic parameters in an anisotropic medium such as a TransverselyIsotropic medium ("TI medium").

BACKGROUND

In geologic exploration and particularly in the exploration forunderground hydrocarbon deposits, it is useful to obtain a measure ofthe elastic parameters of earth strata, either for the direct goal ofunderstanding the behavior of the earth formations subject to mechanicalstresses and strains, or as an intermediate step in improving seismicmeasurement and processing methods.

The elastic moduli of earth formations are directly related to, anduseful for characterizing, rock strength, porosity, lithology, and porefluids.

The elastic moduli of earth strata are also directly related to theacoustic propagation characteristics of the earth, and along withdensity constitute primary parameters of earth strata in seismology. Theelastic moduli are related to the speed of propagation for elastic waveswhich may vary as a function of direction of propagation and directionof polarization. Waves with polarization direction approximately alignedwith propagation direction are referred to as quasi P waves (qP). Waveswith polarization direction approximately orthogonal to propagationdirection are quasi shear waves (qS). The elastic moduli may beexpressed as matrix elements c_(ij) where i≦6, j≦6. The number ofindependent elastic parameters c_(ij) can vary between two parameters inthe case of an isotropic medium and 21 parameters for an arbitraryelastic medium.

The earth is comprised of many layers of geologic strata, and it iscommon to use a horizontally layered medium as a basic model of theearth, wherein each layer is attributed with elastic parameterscorresponding to measured parameters such as the density and thepropagation velocities for P, SV and SH waves. When anisotropy isencountered, it has been found to be particularly useful to model theearth as a TI elastic medium.

A TI medium is transversely isotropic with respect to some symmetryaxis. Using the vertical axis as the symmetry axis, it is referred to asa "TIV medium", whereas if the horizontal axis is used as a symmetryaxis, it is sometimes referred to as a "TIH medium". In either case, aTI medium is characterised by its density ρ together with fiveindependent elastic moduli relating stress and strain. In condensedsubscript notation, with rotational symmetry around the 3-axis, themoduli are {c₁₁,c₁₃,c₃₃,c₅₅,c₆₆ }. The density-normalised moduli, A_(ij)=c_(ij) /ρ, have dimensions of velocity² and are related to elasticpropagation in the medium as follows. Any plane harmonic wavepropagating in the medium, having a phase slowness vector (p₁,p₃) lying,for example, in the 1-3 plane must have displacement either in the 1-3plane (the qP or qSV case) or normal to it (the SH case). In eithercase, the squared phase slowness vector (X,Z)=(p₁ ²,p₃ ²) must satisfythe appropriate Christoffel relation. In the SH case that is

    A.sub.66 X+A.sub.55 Z-1=0;                                 (1)

and in the qP or qSV case that is

    A.sub.11 A.sub.55 X.sup.2 +A.sub.33 A.sub.55 Z.sup.2 +AXZ-(A.sub.11 +A.sub.55)X-(A.sub.11 +A.sub.55)Z+1=0                     (2)

where

    A=A.sub.11 A.sub.33 +A.sub.55.sup.2 -(A.sub.13 +A.sub.55).sup.2 ( 3)

The density normalized moduli in a TI medium correspond to physicalparameters as follows:

    A.sub.11 =A.sub.22 =v.sup.2                                ( 4a)

for horizontally propagating P-waves

    A.sub.33 =v.sup.2                                          ( 4b)

for vertically propagating P-waves

    A.sub.44 =A.sub.55 =v.sup.2                                ( 4c)

for vertically propagating shear waves and horizontally propagatingshear waves with vertical polarization (SV)

    A.sub.66 =v.sup.2                                          ( 4d)

for horizontally propagating shear waves with horizontal polarization(SH)

    A.sub.13 =A.sub.23                                         ( 4e)

is coupled to off-axis qP and qS propagation.

The 3 axis is designated the vertical axis of this TIV medium. Inpractice, where the earth strata of interest exhibit primarilyhorizontal layering, the TIV medium provides a good model for seismicpurposes. If additional parallel vertical fractures are imposed inaddition to the layering, another form called a Fractured TIV or "FTIVmedium" is created, which has especially relevant and useful parallelswith the hydrocarbon bearing earth formations of interest to seismicexploration. A FTIV medium has eight independent elastic parametersc_(ij). Both TI and FTIV media fall into the broader category oforthorhombic media, which are defined as having three mutuallyorthogonal planes of mirror symmetry. A general orthorhombic medium hasnine independent elastic parameters c_(ij).

In laboratory studies the acoustic slowness measurements are sometimesmade in three directions: two axial directions plus one at 45 degrees.See Jones, L. E. A., and Wang, H. F., 1981, "Ultrasonic velocities inCretaceous shales from the Williston basin" Geophysics, 46, 288-297. Insuch a controlled experiment, the axial measurements can be taken todetermine directly A₁₁, A₃₃, A₅₅ and A₆₆ using equations (4). Given thesingle off-axis measurement point XZ, and the predetermined values ofA₁₁, A₃₃, A₅₅, equations (2) and (3) would additionally yield a simpleexpression for A₁₃.

Where measurements are taken in multiple directions, as is often thecase in field studies, and possibly without axial measurements, therehas not been an accepted exact method for combining the multiplemeasurements to estimate elastic parameters except for the SH case.

In attempting to estimate elastic moduli from multiple measurements, onemay solve for the moduli to best fit a given set of n measured phasepoints {(p₁,p₃)_(i) :i=1,n},n>3, of known wave type. This type of estatehas been done for the case of SH data in White, et al., (1983). Equation(1) becomes an over-determined linear system

    A.sub.66 X+A.sub.55 Z=1,

    X= p.sub.1i.sup.2 :i=1,n!,Z= p.sub.3i.sup.2 :i=1,n!,1= 1:i=1,n!(5)

that can be solved for A₆₆ and A₅₅ by standard methods.

The qP and qS cases have commonly been treated by introducingapproximate solutions to equation (2) based on the assumption that amedium is only weakly anisotropic. For example in Gaiser, J. E.,"Transversely isotropic phase velocity analysis from slownessestimates", J. Geophys. Res., 95, 11241-11254 (1990), it has beenproposed to approximate the slowness sureface for qP waves, definedimplicitly by equation (2), by an approximate function having the form

    v.sup.2 ≈C.sub.0 +C.sub.2 cos(2θ)+C.sub.4 cos(4θ)+C.sub.6 cos(6θ)

where in the notation of Equation (2) hereinabove, ##EQU1## There hasbeen no known way of directly solving the Christoffel relation (2)exactly without use of approximations or weak anisotropy assumptions.Also see White, J. E., Martineau-Nicoletis, L., and Monash, C.,"Measured anisotropy in Pierre shale" Geophysical Prospecting,31,709-725 (1983).

Accordingly, it is an object of the invention to determine the elasticmoduli of anisotropic earth strata using acoustic and density datawithout making simplifying approximations of the Christoffel relation orthe assumption of weak anisotropy.

It is also an object of the invention to improve the estimation ofelastic moduli of anisotropic earth formations based on data from wideaperture VSP and density logs.

It is also an object of the invention to improve the speed, accuracy andmethodology in the determination of elastic parameters of earth media.

SUMMARY OF THE INVENTION

In the present invention the elastic moduli are not estimated usingapproximated solutions to the Christoffel relation. Instead, theChristoffel relation is converted algebraically into a linear system andsolved analytically to generate a solution for estimating elasticmoduli. In a preferred embodiment of the invention, a set of qP phasepoints obtained from walkaway VSP data, with multiple vertical angles,plus at least one vertical SV point, are used to derive the elasticmoduli A₁₁,A₁₃,A₃₃,A₅₅. In the case of FTIV or other orthothrombicmedia, the present invention may also be applied to obtain additionalelastic moduli.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is a flow diagram showing the steps in a preferred embodiment ofthe method of the invention in a TI medium.

FIG. 2 shows VSP data used as inputs in accordance with the presentinvention.

FIG. 3 shows vertical apparent slownesses for common source gatherscomputed from the VSP data of FIG. 2.

FIG. 4 shows the horizontal apparent slowness related to the VSP data ofFIG. 2.

FIG. 5 shows a cross plot of estimated qP phase slowness points relatedto FIGS. 2-4, together with the analytic phase slowness curve for thebest fitting TIV medium.

FIG. 6 is a flow diagram showing an extension of the present inventionto determining elastic parameters in a Fractured TIV medium.

DETAILED DESCRIPTION

TI Algebraic Transformation

The qP-qSV Christoffel relation (2) can be rewritten as

    A.sub.11  A.sub.55 X.sup.2 -X!+A.sub.33  A.sub.55 Z.sup.2 -Z!+A XZ!= A.sub.55 (X+Z)-1!.                             (6)

Given (squared) phase data points {(X_(i),Z_(i)):i=1,n} where n=numberof phase points data, and a value for A₅₅, each of the bracketedquantities becomes a data vector and Equation (6) becomes a linearsystem similar to equation (3). With this altered view of theChristoffel relation, it becomes possible to consider other approachesto deriving values for the elastic parameters. Specifically, transformthe original data points (X_(i),Z_(i)) to new variables U_(i), V_(i),W_(i), D_(i)) defined as

    U.sub.i =A.sub.55 X.sub.i.sup.2 -X.sub.i, V.sub.i =A.sub.55 Z.sub.i.sup.2 -Z.sub.i, W.sub.i =X.sub.i,Z.sub.i, D.sub.i =A.sub.55 (X.sub.i +Z.sub.i)-1.                                              (7)

Let

U= U_(i) : i=1, . . . ,n!

V= V_(i) : i=1, . . . ,n!

W= W_(i) : i=1, . . . ,n!

D= D_(i) : i=1, . . . ,n!

Then equation (6) can be rewritten

    A.sub.11 U+A.sub.33 V+AW=D                                 (8)

Equation (8) represents a linear system having the unknown scalarcoefficients A₁₁, A₃₃, and A, and is exactly equivalent to the statementthat each data point satisfies the Christoffel relation. In matrixnotation this takes the form ##EQU2##

A solution to such a linear problem can be computed exactly and quickly,for example using a commercial package such as the Mathematica Programavailable from Wolfram Research. Inc. Cf. Mathematica (Addison Wesley,1971) pp. 659-662. Other representations and equivalent solutions to theabove linear system will be obvious to skilled persons in this field,and will not be discussed further. Once A₁₁, A₃₃ and A have beendetermined, the elastic parameter A₁₃ can then be obtained from equation(3) assuming A₁₃ +A₅₅ >0, using the following equivalent form:

    A.sub.13 =(A.sub.11 A.sub.33 +A.sub.55.sup.2 -A).sup.5 -A.sub.55. (10)

An alternative solution for A₁₃ (with A₁₃ +A₅₅ <0) would yield a mediumwith the same phase slowness surfaces, but with anomalous polarizationsnear 45 degrees. For inversion of physical measurements, this case isunlikely to occur and can easily be distinguished by looking atpolarizations near 45 degrees. Mathematically, it is simply a secondsolution to the inversion problem which is physically permitted providedcertain stability conditions are met.

With the above understanding of the reparametrization, and referring tothe figures, a preferred form of the method will be described.

Referring to FIG. 1, a value for the parameter A₅₅ is estimated orcomputed, at step 15, using any of various known methods. For example itis known to derive A₅₅ from SH Analysis 13 (See White, 1983), or fromFree Parameter Analysis 14 indexing of a family of solutions to beoptimised by a one-parameter search. It is also possible to derive A₅₅from the Sonic Log 11 as the velocity-squared of the verticallypropagating shear wave, or from Vertical Seismic Profile (VSP) data 12as described below.

Raw VSP data 12 are shown in FIG. 2 as a function of seismic sourceoffset (in kilometers). The vertical apparent slownesses are calculatedfor each common source gather and shown in FIG. 3 also as a function ofsource offset. The data points near 0 offset are missing because thedrilling rig, located at point 0--0, prevented seismic shots from beingtaken at that point. The four vertical slowness curves are, from top tobottom, the downgoing qS wave 22, downgoing qP wave 23, reflected qPwave 24, and reflected qS wave 25, respectively. The 0 offset interceptpoint 26 can be used in Equation (4c) to obtain an estimate of theparameter A₅₅ at step 15, as mentioned above. The downgoing qP wavecurve 23 is the requisite vertical slowness ∂T/∂z as a function ofsource offset.

The VSP data 12 are also used to generate the horizontal apparentslowness 30 shown in FIG. 4 as a function of offset. Referring to FIG. 1step 17, the VSP data 12 is used to estimate the travel time T(x) andhence ∂T/∂x, the apparent qP horizontal slowness. The vertical slowness∂T/∂z, and horizontal slowness ∂T/∂x, each as a function of offset, arecombined to generate the cross plot of phase slowness points curve 40shown in FIG. 5.

The Christoffel relation is reparametrized at 19 and the phase slownesspoints 40 are used as input data to transform to the Linear System 20defined in equation (8). The phase slowness data 40 comprise a largenumber of phase slowness points scattered over almost the entire rangeof vertical angles. In fact, in accordance with the present invention,it is not necessary to obtain such a large amount of high quality phaseslowness points, and only a few, for example three points covering awide range of angles, would be adequate. This yields the elasticparameters A₁₁, A₁₃, and A₅₅ at step 21, as explained above. Curve 41shown in FIG. 5 is the analytic slowness curve associated with theelastic moduli A₁₁, A₁₃, A₃₃, and A₅₅ that were obtained using themethod herein described.

The estimation of the last independent elastic parameter for a TIVmedium, A₆₆, is already known, for example using the technique describedin Ellefsen, K. J. et al, "Estimating a shear modulus of a transverselyisotropic formation", Geophysics, Vol. 57, No. 11 (November 1992) p.1428-1434, or White (1983) mentioned earlier.

The method described above can be extended to use other types ofacoustic phase point data, and to generate elastic parameters in mediaother than the TI type.

In a horizontally stratified offshore environment, offset VSP data canbe acquired along multiple azimuthal lines to obtain, for each measuredazimuth, qP phase points at all vertical angles, together with a singlevertical SV point. With this type of measurement set, and consideringthat the true medium may show azimuthal as well as vertical anisotropy,it is possible to use the TI Transformation in a more general mediumwhich is not simply Transversely Isotropic. For example, given varioussymmetry conditions of a non-TI medium, it will be possible to use themethod of the present invention to obtain elastic parameters, whiledeviating from the detailed preferred embodiment shown and describedherein.

Given data from a single azimuth, as in the example discussed earlier,or from multiple azimuths in which the data is seen to be independent ofazimuth, density-normalised moduli characterising the medium may bedirectly obtained as described. When the data show variation withazimuth, provided certain symmetry conditions are met, the TITransformation can be used as part of a method to obtain the elasticmoduli for the medium. For example, the determination of elasticparameters in a Fractured TIV medium using the present invention isshown in FIG. 6, and will be described below.

A FTIV medium, assuming a single fracture set aligned with the 2-axishas elastic moduli that satisfy the symmetry relations of anorthothrombic medium, together with an additional relation given inHood, J. A., and Schoenberg, M., "Estimation of vertical fracturing frommeasured elastic moduli,"Geophys. Res., 94, 15611-15618 (equation 13)(1989):

    A.sub.12 =(A.sub.13 A.sub.22 -A.sub.11 A.sub.23)/(A.sub.23 -A.sub.13). (11)

It is known that for qP and qSV propagation with the slowness vector inone of the symmetry planes, the Christoffel relation reduces to the TIform of Equation (2). In the 1-3 plane, Equation (2) holds as written.In the 2-3 plane, in Equation (2), one must change all subscript 1's to2's and 5's to 4's. In the 1-2 plane, in Equation (2), one must changeall subscript 3's to 2's and 5's to 6's. Based on this transformationand the above observation, the inversion problem for FTIV media isreduced to the inversion problem for these associated TI cases.

The data set that is needed, assuming knowledge of the symmetrydirections, would include phase slowness measurements, for exampleWalkaway VSP ("WVSP"), in the two symmetry planes plus a third at anintermediate azimuth, for example at about 45 degrees thereto. Apractical sequence for obtaining the elastic parameters of a FracturedTIV medium would comprise the following.

First, identify the symmetry directions and determine A₅₅ and A₄₄. Thismay be done in advance using shear sonic measurements, or by examiningmultiple azimuth walkaway VSP data. An estimate for vertical qSVslowness can be obtained in a single azimuth by observing upgoing shearwaves at near vertical incidence. The 2-3 plane (parallel to thefractures) can be determined from multi-azimuth data as the azimuth withthe minimum qP slowness at all vertical angles. This information couldalso be obtained from polarization analysis of vertical shear waves(shear splitting) in shear VSP or shear sonic data. If the data isindependent of azimuth, then the medium is TIV and can be handled by theuse of the TI Algebraic Transformation and Linear System solution aspreviously described. If the medium does not show two orthogonal planesof mirror symmetry, then this is not an orthorhombic medium and it isnoted that the present method would not apply exactly.

Next, use the TI Algebraic Transformation to invert the qP data in the1-3 plane for {A₁₁, A₁₃, A₃₃, A₅₅ } and in the 2-3 plane for{A₂₂,A₂₃,A₃₃,A₄₄ }. The elastic parameter A₁₂ may be derived directlyusing Equation (11).

Finally, invert a set of at least three qP points from the 1-2 plane todetermine A₆₆. This can be done by using the TIV transformation appliedto the 1-2 plane qP points to determine A₁₂, and altering the parameterA₆₆ in the transformation until the output value of A₁₂ satisfiesEquation (11).

In a more general orthorhombic medium, Equation (11) may not hold.However, in situations where an independent estimate of A₆₆ can be made(say, from use of tube-wave data, or in a lab study), then thehorizontal qP points can be used to determine A₁₂ using the TITransformation in the 1-2 plane as in Step (b) without use of Equation(11). If only vertical shear points are available and the medium is notfar from a fractured TIV medium (in the sense that Equation (11) gives aplausible value for A₁₂), then the method described in the previousparagraph is likely to yield a fractured TIV medium that agrees toexperimental accuracy with the true medium at all the measured phasepoints.

The elastic parameters determined in the present invention are usefulfor various aspects of seismic processing. The elastic moduli contributeto the velocity models used for anisotropic migration, and serve tocalibrate AVO studies. Elastic moduli are also fundamental properties ofrocks that are important in determining and controlling wellborestability, permeability, shale maturation, pore pressure and otherrelated parameters.

It is clear that the methods described herein are not limited to use ofthe particular mathematical forms and conventions used, and that theelastic parameters Aij are essentially equivalent to the parameters Cijused in the cited publications, and these are also convertible to otherrepresentations. Similarly, it is obvious that the mathematicalexpressions have various equivalents, which may be encompassed withinthe scope of the present invention.

I claim:
 1. A method of characterizing anisotropic properties of anearth medium, comprising the steps of:a) acquiring phase slowness datawhich satisfy a Christoffel relationship from seismic measurements ofthe earth medium; b) transforming the phase slowness data so as toprovide a series of equivalent linear equations which linearize theChristoffel relationship; c) deriving a solution of the system of linearequations resulting from the transformation; and d) deriving the elasticproperties of the earth medium from the solution.
 2. A method as claimedin claim 1, wherein the step of transforming the phase slowness dataprovides a series of linear equations which linearize a mathematicallyequivalent representation of the Christoffel relationship.
 3. A methodas claimed in claim 2, wherein the phase slowness data are transformedaccording to

    U= A.sub.55 X.sub.i.sup.2 -X.sub.i :i=1,n!, V= A.sub.55 Z.sub.i.sup.2 -Z.sub.i : i=1, n!, W= X.sub.i Z.sub.i :i=1,n!, D= A.sub.55 (X.sub.i +Z.sub.i)-1:i=1, n!

and using the system of linear equations defined by

    A.sub.11 U+A.sub.33 V+AW=D

to represent the Christoffel relationship from which the earth elasticproperties are derived.
 4. A method as claimed in claim 1, wherein atleast three phase slowness points are used for transforming the data soas to characterize the anisotropic properties of a TI earth medium.
 5. Amethod as claimed in claim 1, comprising using an independent estimationof elastic parameter A₅₅ as an input to the transforming step.
 6. Amethod as claimed in claim 3, comprising using an independent estimationof elastic parameter A₅₅ as an input to the transforming step.
 7. Amethod as claimed in claim 3, further comprising obtaining phaseslowness data from at least one symmetry plane of an orthorhombic mediumand applying the transformation to the phase slowness data to deriveelastic parameters associated with the at lest one symmetry plane.
 8. Amethod as claimed in claim 7, comprising varying A₆₆ as a test parameterto generate a value of A₁₂ in accordance with the transforming step tomatch a value of A₁₂ generated by A₁₂ =(A₁₃ A₂₂ -A₁₁ A₂₃)/(A₂₃ -A₁₃). 9.A method as claimed in claim 1, further comprising generating a velocitymodel of the earth medium using the elastic parameters.
 10. A method asclaimed in claim 1, further comprising generating an AVO model of theearth medium using the elastic parameters.
 11. A method as claimed inclaim 1, further comprising generating a model for mechanical propertiesof the earth medium using the elastic parameters.